///////////////////////////////////////////////////////////////////////////////////////////////////
/// \file correlator.h
///
/// \brief   header file for the correlator class
///
///////////////////////////////////////////////////////////////////////////////////////////////////

#ifndef PCPS_CORRELATOR_HEADER
#define PCPS_CORRELATOR_HEADER

#include<src/pcps.h>
#include<src/input.h>
#include<src/coeff_obj.h>
#include<src/mantissa_exponent.h>
#include<src/wfn_data.h>
#include<src/command.h>

namespace pcps {

  ///////////////////////////////////////////////////////////////////////////////////////////////////
  /// \brief   Class to represent a correlator.
  ///
  ///////////////////////////////////////////////////////////////////////////////////////////////////
  template<class S> class Correlator : public pcps::CoeffObj<S> {

    public:

      /// \brief function to access the index data object by reference
      pcps::WfnData<int> & index_data() { return *this->_wi_deps->at(0); }

      /// \brief function to access the index data object by constant reference
      const pcps::WfnData<int> & index_data() const { return *this->_wi_deps->at(0); }

      /// \brief function to access the value data object by reference
      pcps::WfnData<S> & value_data() { return *this->_wf_deps->at(0); }

      /// \brief function to access the value data object by constant reference
      const pcps::WfnData<S> & value_data() const { return *this->_wf_deps->at(0); }

      /// \brief function to access the indices by reference
      std::vector<int> & indices() { return this->index_data().data_vec(); }

      /// \brief function to access the indices by constant reference
      const std::vector<int> & indices() const { return this->index_data().data_vec(); }

      /// \brief function to access the values by reference
      std::vector<S> & values() { return this->value_data().data_vec(); }

      /// \brief function to access the values by constant reference
      const std::vector<S> & values() const { return this->value_data().data_vec(); }

      /// \brief function to return a pointer to the values
      S * data_ptr() {
        assert( this->values().size() > 0 );
        return &this->values()[0];
      }

      /// \brief function to return a pointer to the values
      const S * data_ptr() const {
        assert( this->values().size() > 0 );
        return &this->values()[0];
      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief Constructs the correlator by reading it from an archive.
      ///
      /// \param[in]      userinp  the user's input options
      /// \param[in,out]  arch     an archive from which to read the correlator
      /// \param[in]      wi_map   map of names to available integer wavefunction data objects
      /// \param[in]      wf_map   map of names to available floating point wavefunction data objects
      /// \param[in]      co_map   map of names to available coefficient objects
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      Correlator(const pcps::Input & userinp,
                 pcps::Archive & arch,
                 const std::map<std::string, boost::shared_ptr<pcps::WfnData<int> > > & wi_map,
                 const std::map<std::string, boost::shared_ptr<pcps::WfnData<S> > > & wf_map,
                 const std::map<std::string, boost::shared_ptr<pcps::CoeffObj<S> > > & co_map)
        : CoeffObj<S>(1, 1, 0)
      {

        // get name of the correlator
        pcps::mpi::read_and_bcast(arch, this->_name, "Failed to read name in Correlator constructor.");

        // get the name of the indices data
        std::string wi_name;
        pcps::mpi::read_and_bcast(arch, wi_name, "Failed to read name of index data in Correlator constructor.");
        if (wi_map.count(wi_name) == 0)
          throw pcps::Exception("no integer data by the name %s") % wi_name;
        this->_wi_deps->at(0) = wi_map.find(wi_name)->second;

        // check that all indices are sane
        if (pcps::mpi::rank() == 0)
          for (std::vector<int>::const_iterator ind = this->indices().begin(); ind != this->indices().end(); ind++)
            if ( *ind < 0 || *ind >= userinp.nsites() )
              throw pcps::Exception("impossible correlator site \"%i\" in Correlator constructor.") % (*ind);

        // get values
        std::string wf_name;
        pcps::mpi::read_and_bcast(arch, wf_name, "Failed to read name of value data in Correlator constructor.");
        if (wf_map.count(wf_name) == 0)
          throw pcps::Exception("no floating point data by the name %s") % wf_name;
        this->_wf_deps->at(0) = wf_map.find(wf_name)->second;

        // check that the value data object has the correct number of elements
        if ( this->values().size() != pcps::pow<size_t>(userinp.configs_per_site(), this->indices().size()) )
          throw pcps::Exception("Wrong number of values in Correlator constructor.  Found %i values, expected %i.")
                % this->values().size() % pcps::pow<size_t>(userinp.configs_per_site(), this->indices().size());

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief converts the object into the command that will create the object
      ///
      /// \param[in,out]  arch     an archive in which to store the command
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void to_command(pcps::Archive & arch) const {
        arch << pcps::CorrelatorCommand::keyword() << std::endl;
        arch << this->_name                        << std::endl;
        arch << this->index_data().name()          << std::endl;
        arch << this->value_data().name()          << std::endl;
        arch << std::endl;
      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief returns the value of the coefficient for an arbitrary configuration
      ///
      /// \param[in]      config   the configuration whose coefficient we want
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      pcps::MantExp<S> get_value(const pcps::Configuration<S> & config) {
        throw pcps::Exception("NOT YET IMPLEMENTED");
        return pcps::MantExp<S>();
      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief returns the value of the coefficient for a configuration that differs from the
      ///        currently initialized configuration as specified.
      ///
      /// \param[in]      new_config  the lattice configuration whose coefficient we want
      /// \param[in]      old_config  the lattice configuration the object is currently initialized for
      /// \param[in]      diff        an object specifying how the configurations differ
      ///
      /// \return the coefficient at the new configuration
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      pcps::MantExp<S> get_value(const pcps::Configuration<S> & new_config,
                                 const pcps::Configuration<S> & old_config,
                                 const pcps::ConfigDiff<S> & diff)
      {
        throw pcps::Exception("NOT YET IMPLEMENTED");
        return pcps::MantExp<S>();
      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief Computes the ratios of the noiseless coefficients at each singly excited configuration
      ///        to the noisy coefficient at the current configuration
      ///
      /// \param[in]      config      the lattice configuration the object is currently initialized for
      /// \param[out]     a_ratios    on exit, the matrix of excitation ratios for alpha excitations
      /// \param[out]     b_ratios    on exit, the matrix of excitation ratios for beta excitations
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void single_excitation_ratios(const pcps::Configuration<S> & config,
                                    std::vector<S> & a_ratios,
                                    std::vector<S> & b_ratios)
      {
        throw pcps::Exception("NOT YET IMPLEMENTED");
      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief Adds the derivatives of the coefficient with respect to the wavefunction's data
      ///        variables to the supplied vector after multiplying the derivatives by a constant.
      ///        The derivatives are computed assuming that the object is currently initialized for the
      ///        supplied lattice configuration.
      ///
      /// \param[in]      config  the lattice configuration the object is currently initialized for
      /// \param[in]      alpha   constant to multiply the coefficient's derivatives by
      /// \param[out]     der_vec vector to which to add the derivative contributions from this object
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void add_derivatives(const pcps::Configuration<S> & config,
                           const pcps::MantExp<S> & alpha,
                           std::vector<S> & der_vec)
      {
        throw pcps::Exception("NOT YET IMPLEMENTED");
      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief For each orbital pair excitation, contracts the wavefunction's derivatives at the
      ///        corresponding excited configuration with the supplied trial vector.
      ///        It is assumed that the object is initialized for the supplied configuration.
      ///
      /// \param[in]      config      the lattice configuration to excite from.
      /// \param[in]      alpha       constant to multiply the contraction by
      /// \param[in]      trial_vec   the vector to be contracted
      /// \param[in,out]  out_vec_a   vector to which to add the alpha orbital result
      /// \param[in,out]  out_vec_b   vector to which to add the beta orbital result
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void contract_der_to_orb_pairs(const pcps::Configuration<S> & config,
                                     const pcps::MantExp<S> & alpha,
                                     const std::vector<S> & trial_vec,
                                     std::vector<S> & out_vec_a,
                                     std::vector<S> & out_vec_b)
      {
        throw pcps::Exception("NOT YET IMPLEMENTED");
      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief For each wavefunction variable, contracts the various singly excited configurations'
      ///        wavefunction derivatives with respect to that variable with the supplied trial vector.
      ///        It is assumed that the object is initialized for the supplied configuration.
      ///
      /// \param[in]      config        the lattice configuration to excite from.
      /// \param[in]      alpha         constant to multiply the contraction by
      /// \param[in]      trial_vec_a   the vector to be contracted for alpha orbitals
      /// \param[in]      trial_vec_b   the vector to be contracted for beta orbitals
      /// \param[in,out]  out_vec       vector to which to add the result of the contraction
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void contract_orb_pairs_to_der(const pcps::Configuration<S> & config,
                                     const pcps::MantExp<S> & alpha,
                                     const std::vector<S> & trial_vec_a,
                                     const std::vector<S> & trial_vec_b,
                                     std::vector<S> & out_vec)
      {
        throw pcps::Exception("NOT YET IMPLEMENTED");
      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief Multiply the supplied trial vector by the wavefunction derivative outer product matrix
      ///        (the local contribution to the overlap matrix) for the given configuration.
      ///        It is assumed that the object is initialized for the supplied configuration.
      ///
      /// \param[in]      config      the lattice configuration to excite from.
      /// \param[in]      alpha       constant to multiply the contraction by
      /// \param[in]      trial_vec   the vector to be contracted
      /// \param[in,out]  out_vec     vector to which to add the result of the contraction
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void contract_der_to_der(const pcps::Configuration<S> & config,
                               const pcps::MantExp<S> & alpha,
                               const std::vector<S> & trial_vec,
                               std::vector<S> & out_vec)
      {
        throw pcps::Exception("NOT YET IMPLEMENTED");
      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief Initializes the object for the given lattice configuration.
      ///
      /// \param[in]      config  the lattice configuration we want the object initialized for
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void initialize(const pcps::Configuration<S> & config) {
        throw pcps::Exception("NOT YET IMPLEMENTED");
      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief Updates the object so that it is initialized for a new configuration that differs from
      ///        the currently initialized configuration as specified.
      ///
      /// \param[in]      new_config  the lattice configuration we want to object to be initialized for
      /// \param[in]      old_config  the lattice configuration the object is currently initialized for
      /// \param[in]      diff        an object specifying how the configurations differ
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void update(const pcps::Configuration<S> & new_config,
                  const pcps::Configuration<S> & old_config,
                  const pcps::ConfigDiff<S> & diff)
      {
        throw pcps::Exception("NOT YET IMPLEMENTED");
      }

  };

}

#endif
